#ifndef RSIM_SVD_H
#define RSIM_SVD_H

/* -------------------------------------------------------------------------- *
 * File: Svd.h                                                                *
 * Authors: Kishor Bhalerao                                                   *
 * Email : kishor8dm@gmail.com                                                *
 * Contributors:                                                              *
 * Email:                                                                     *
 * copyright (c) 2010 Authors.                                                *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify it    *
 * under the terms of the GNU General Public License as published by the Free *
 * Software Foundation, either version 3 of the License, or any later version.*
 *                                                                            *
 * This program is distributed in the hope that it will be useful, but        *
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
 * for more details. (http://www.gnu.org/licenses/)                           *
 * -------------------------------------------------------------------------- */

#include "Matrix.h"
#include "Vector.h"
#include "LinAlgDecl.h"

#include <algorithm>
#include <cmath>

using namespace std;

namespace RSIM
{
/** Singular Value Decomposition.

<P>
Note: This file was obtained from Template Numerical Toolkit (TNT).
Details of TNT can be found on the website http://math.nist.gov/tnt/
This file is minimally modified to make it compatible with RSIM.
The originial file name (available for download on TNT website) 
is jama_svd.h.

<P>
For an m-by-n matrix A with m >= n, the singular value decomposition is
an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
an n-by-n orthogonal matrix V so that A = U*S*V'.
<P>
The singular values, sigma[k] = S[k][k], are ordered so that
sigma[0] >= sigma[1] >= ... >= sigma[n-1].
<P>
The singular value decompostion always exists, so the constructor will
never fail.  The matrix condition number and the effective numerical
rank can be computed from this decomposition.

<p>
        (Adapted from JAMA, a Java Matrix Library, developed by jointly 
        by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
*/
class SVD{
        private:
                Matrix U, V;
                Vector s;
                int m, n;
                bool isFlat;
        public:
                SVD(const Matrix& Arg);

                void getU(Matrix& A);

                /** Return the right singular vectors */
                void getV(Matrix &A);

                /** Return the one-dimensional array of singular values */
                void getSingularValues(Vector &x){x = s;}
                
                void getColumnBasis(Matrix& M);

                /** Return the diagonal matrix of singular values
                @return     S
                */
                void getS(Matrix &A);

                /** Two norm  (max(S)) */
                double norm2(){return s[0];}

                /** Two norm of condition number (max(S)/min(S)) */
                double cond(){return s[0]/s[min(m,n)-1];}

                /** Effective numerical matrix rank
                @return     Number of nonnegligible singular values.
                */
                int rank();
};

} // namespace RSIM
#endif
